
namespace Tuna {

template<class Tprec, int Dim>
inline 
bool CDS_LES<Tprec, Dim>::calcCoefficients(const ScalarField &nut) { 
    Tprec dyz = dy * dz, dxz = dx * dz, dxy = dx * dy;
    Tprec dyz_dx = dyz / dx, dxz_dy = dxz / dy, dxy_dz = dxy / dz;
    Tprec ce, cw, cn, cs, cf, cb;
    Tprec nutinter;
    Tprec dxyz_dt = dx * dy * dz / dt;

    Tprec Prt = Gamma;

    for (int i =  bi; i <= ei; ++i)
	for (int j = bj; j <= ej; ++j)
	    for (int k = bk; k <= ek; ++k)
	    {
		ce = u(i  , j, k) * dyz;
		cw = u(i-1, j, k) * dyz;		
		cn = v(i, j  , k) * dxz;
		cs = v(i, j-1, k) * dxz;		
		cf = w(i, j, k) * dxy;
		cb = w(i, j, k-1) * dxy;		

/*********************************************************************
// Actually, nut is calculated on center of volumes!!!
//		
// kappa_turbulent = nu_turbulent / Pr_turbulent.
// Because nut is calculated on nodes, we have to interpolate to the
// center of volume.
//
//		nutinter = 0.125 * ( nut(i  ,j  ,k  ) + nut(i  ,j  ,k-1) +
//				     nut(i  ,j-1,k-1) + nut(i  ,j-1,k  ) +
//				     nut(i-1,j  ,k  ) + nut(i-1,j  ,k-1) +
//				     nut(i-1,j-1,k-1) + nut(i-1,j-1,k  ) );


		aE (i,j,k) = (1 + nutinter / Prt) * dyz_dx - ce * 0.5;
		aW (i,j,k) = (1 + nutinter / Prt) * dyz_dx + cw * 0.5;
		aN (i,j,k) = (1 + nutinter / Prt) * dxz_dy - cn * 0.5;
		aS (i,j,k) = (1 + nutinter / Prt) * dxz_dy + cs * 0.5;
		aF (i,j,k) = (1 + nutinter / Prt) * dxy_dz - cf * 0.5;
		aB (i,j,k) = (1 + nutinter / Prt) * dxy_dz + cb * 0.5;
*********************************************************************/

		aE (i,j,k) = (1 + nut(i,j,k) / Prt) * dyz_dx - ce * 0.5;
		aW (i,j,k) = (1 + nut(i,j,k) / Prt) * dyz_dx + cw * 0.5;
		aN (i,j,k) = (1 + nut(i,j,k) / Prt) * dxz_dy - cn * 0.5;
		aS (i,j,k) = (1 + nut(i,j,k) / Prt) * dxz_dy + cs * 0.5;
		aF (i,j,k) = (1 + nut(i,j,k) / Prt) * dxy_dz - cf * 0.5;
		aB (i,j,k) = (1 + nut(i,j,k) / Prt) * dxy_dz + cb * 0.5;
		aP (i,j,k) = aE (i,j,k) + aW (i,j,k) + 
		             aN (i,j,k) + aS (i,j,k) +
		             aF (i,j,k) + aB (i,j,k) + dxyz_dt;
//		+ (ce - cw) + (cn - cs) + (cf - cb);

		sp (i,j,k) = phi_0(i,j,k) * dxyz_dt ;
	    }

    applyBoundaryConditions3D();
    return 1;
}

} // Tuna namespace















